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Abstract. We investigate the Bayesian framework for detection of continuous 
gravitational waves (GWs) in the context of targeted searches, where the phase 
evolution of the GW signal is assumed to be known, while the four amplitude 
parameters are unknown. We show that the orthodox maximum-likelihood 
statistic (known as J^-statistic) can be rediscovered as a Bayes factor with an 
unphysical prior in amplitude parameter space. We introduce an alternative 
detection statistic ( "S-statistic" ) using the Bayes factor with a more natural 
amplitude prior, namely an isotropic probability distribution for the orientation of 
GW sources. Monte-Carlo simulations of targeted searches show that the resulting 
Bayesian £>-statistic is more powerful in the Neyman-Pearson sense (i.e. has a 
higher expected detection probability at equal false-alarm probability) than the 
frequentist jF-statistic. 



PACS numbers: 02.50.Tt,02.70.Rr,04.30.w,07.05.Kf,95.85.Sz 
1. Introduction 

Searches for gravitational waves (GWs) often consist of testing the data for the 
presence of signals from a known family of waveforms, parametrized by (generally 
unknown) signal parameters. Here we consider the class of coherent GW signals of 
constant amplitude and polarization, which include "continuous GWs", e.g. from non- 
axisymmetric spinning neutron stars (see [1] for a review) , stellar-mass binary systems 
in the LISA frequency band (e.g. [2J|3]), and coalescence of (non-precessing) binary 
systems [J. 

We distinguish two classes of signal parameters: (i) four "amplitude parameters", 
namely the amplitudes A + and A x of the two GW polarizations, the orientation angle 
tp of the principal polarization axis, and the initial GW phase (f>o, and (ii) the set of 
"Doppler parameters" A, which determine the time evolution of the GW phase cf>(t; A). 
We restrict our attention to targeted searches, in which the Doppler parameters A are 
assumed to be known, resulting in a detection problem with four unknown amplitude 
parameters. 

The popular, yet ad-hoc, orthodox approach consists of maximizing the likelihood 
function of the data over these four amplitude parameters. It was first shown in [5] 
that this maximization can be achieved analytically, resulting in a computationally 
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very efficient detection statistic, known as the ^"-statistic, which has been used in a 
number of searches for GWs (e.g. [51 [7J [HI IS]). 

Here we investigate an alternative Bayesian approach, which leads us to the Bayes 
factor as a useful classical detection statistic (see also [TUl HH HU H3] )• Contrary 
to the maximum-likelihood approach, the Bayesian framework requires the explicit 
prescription of a prior probability distribution for the unknown signal parameters. 
We show that a particularly simple, yet unphysical, choice of amplitude prior results 
in the JF-statistic as a special case of a Bayes factor. This illustrates that frequentist 
ad-hoc statistics often carry their own unchecked and implicit priors, hidden from view 
and often unknown to the user (see also O [TS] ) . 

We can derive a more natural amplitude prior from our model assumption 
about the emission of GWs from non-axisymmetric spinning systems: the amplitude 
parameters are closely related to the orientation of the emitter with respect to 
the observer. In the absence of astrophysical information, an isotropic probability 
distribution for the spin-axis orientation is therefore the natural choice. We refer to 
the Bayes factor resulting from this amplitude prior as the "S-statistic" . Isotropic spin- 
axis orientation priors have been used previously for Bayesian parameter estimation 
[51[TH1[T7], and in Monte-Carlo simulations to determine frequentist upper limits [71 ITS], 
A comparison of Bayesian and frequentist methods for setting upper limits can be 
found in PH©. 

The JF-statistic (or matched filtering in general) has often been incorrectly referred 
to as an "optimal statistic" (e.g. [5j [6j [H [20j [9] to name only a few). Using Monte- 
Carlo simulations we show that the Bayesian ^-statistic is more powerful (i.e. has a 
higher expected detection probability at equal false-alarm probability) than the J-- 
statistic for GWs emitted by systems with random (isotropic) spin-axis orientations. 
This is a direct consequence of the B-statistic prior being consistent with the injected 
distribution of parameters, contrary to the (implicit) .F-statistic prior. Similar results 
were found previously in the case of burst detection statistics [12j [ill US] ■ 



2. Signal model: coherent gravitational waves of constant amplitudes 

The spatial metric perturbation h(t) of a coherent GW of constant amplitudes (far 
from the source) can be written as 

h(t) = e+ A+ cos [4>(t; A) + <j> ] + e x A x sin [cf>(t; A) + fa] , (1) 

where *e + = u ® u — v ® v, and *e x = u<g> v + v (3 ii are two polarization basis tensors, 
constructed from a right-handed basis {«, v, — h}. The unit-vector h is pointing along 
the line of sight from the detector to the source, and the wave-plane basis vectors 
{u, i>} are aligned with the principal polarization axes of the GW. In general the GW 
phase (f>(t; A) depends on a set of Doppler parameters A, which include the source 
sky-position n, the GW frequency /, and possibly higher-order time derivatives of the 
frequency {/, /, . . .}. If the source is a neutron star in a binary system, then A would 
also include the orbital parameters of the system. 

We assume the GW emitter consists of a (non-axisymmetric) rotating system 
with spin f2 and ellipticity e with respect to the rotation axis. The corresponding 
characteristic amplitude ho of the GW at the detector can be expressed as 

, 4G/ zz ft 2 

ho = ^^s, (2) 
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where I zz is the moment of inertia with respect to the rotation axis, the rotation 
rate is f2 = |f2|, and d is the distance to the detector. The emitter geometry is 
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Figure 1. Source geometry angles: i is the inclination angle between the rotation 
axis Q and the line of sight (—ft)- The polarization angle tp measures the 
orientation of the projected rotation axis Qy in the sky plane with respect to 
an observer frame {§, £}. 

fully characterized by two Euler angles describing the orientation of the rotation axis, 
namely the inclination angle i of the rotation axis with respect to the line of sight h, 
and the polarization angle ip describing the orientation of the projected rotation axis 
in the sky-plane (see Fig. [I]). The two polarization amplitudes A + , A x in can be 
expressed as 

A + = ^o(l + cos 2 t) , and A x =IiqCosl, (3) 

in terms of the characteristic amplitude hg and the inclination angle l. This 
corresponds to a choice of the wave-plane coordinate axis u perpendicular to the 
rotation axis [21 and assigns A + to the larger of the two principal polarization 
amplitudes, i.e. A + > A x . For each sky position n we can define a source-independent 
(right-handed, orthonormal) basis {£, Cj — n}, e.g. where £ lies in the ecliptic plane and 
£ in the northern hemisphere. This allows us to define the polarization angle ip as the 
angle between the principal polarization axis u of the GW and the basis vector £, i.e. 
ip = Z(u,£) = Z(v,Q, where it oc f2y (see Fig. lk. As first shown in [5], the strain h(t) 
measured in the detector due to a GW signal iQ can be expressed as 

h(t;A,\)=A»h fl (t;\), (4) 

where we use automatic summation over repeated amplitude indices fj,. The 

explicit form of the four "basis functions" h^(t;X) can be found in [9], for example, 
but is not important for the following discussion. The four canonical amplitudes A 11 
are defined as 

A 1 = A + cos 2^ cos 0o ^ A x sin 2-0 sin^o , 

A 2 = A + sin 2-0 cos <po + A x cos 2-0 sin <p$ , 

, (5) 
A = —A + cos 2tp sin (f> — A x sin 2ip cos tp , 

A 4 = —A + sin.20 sin0 o + A x cos 2ip cos^o . 
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The set of amplitude parameters A can therefore be expressed cither in "physical 
coordinates" A 1 , i.e. 

{A} 1 = A l = {/i ,cost,V,<M, (6) 
or in "canonical coordinates" A 11 given in (p>]). 



3. Hypothesis testing for GW detection 

3.1. Simple versus composite hypotheses 

In the following we restrict ourselves to the case where we need to decide only between 
two hypotheses, namely Ti^ = "the data x consist of only noise n" , and Tis = "the 
data contains a signal s in addition to noise n" , where we assume a signal s — h(t; A, A) 
of the form Q . For simplicity we assume the Doppler parameters A = A s to be known 
a-priori, while the four amplitude parameters A € A are unknown. This corresponds 
to a targeted search for GWs, for example from an isolated pulsar with known sky 
position and GW frequency f(t). We can formally write the two hypothesis as 

TC N : x(t) = n{t) , 

Tis : x(t) = n(t) + s(t; A, A s ) , for any A G A . 
Note that 7Yn is a simple hypothesis, which means that all model parameters are fully 
specified, namely s — 0. On the other hand, Tis is a composite hypothesis, as the 
amplitude parameters A G A of the signal are unknown. The composite hypothesis 
TLs can be considered as a union of simple hypotheses, i.e. TLs = U TLs(A). 



3.2. Scalar product and Gaussian noise 

Assuming Gaussian stationary noise with known (single-sided) power spectral density 
5 n , the probability density function (pdf) for a particular pure-noise time series x = n 
can be written as 

pdf (n|5 n ) - «^ ( " |n) , (8) 

where k is a normalization constant and we defined the scalar product (x\y) between 
time-series x(t) and y(t) as 

(x\y)=4Rj" i{ sJ f \ f) df, (9) 

where x(f) denotes the Fourier transform of x(t) and * denotes complex conjugation 
(e.g. see [H]). The likelihood of observing data x(t) in the presence of a signal s(t) is 
therefore 

pdf {x\sS a ) = Ke -i« x - s ^ x - s » . (10) 



3.3. Neyman- Pearson optimality 

A detection statistic &(%) is a real- valued function of the data x, such that the 
corresponding test of threshold 6* decides for Hn if 6(x) < 6*, and for Tts if 
&(x) > 6*. Such a test is typically subject to two types of error: a "false alarm", i.e. 
choosing TCs when TCn is true, and a "false dismissal", i.e. rejecting Tts when it is in 
fact true. We can express the probability of a false alarm as 

poo 

f A (6*) = P(6>6*\H N )= / P df(6|W N )d6. (11) 

Je, 
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The probability fu of a false dismissal is fn(&*; IIa) — P(& < S*|Hs IIa), where we 
notice that the specification ^ of the signal hypothesis Tis is in fact incomplete: in 
addition to the condition A € A, we also need to specify the probability density for A, 
which we denote as IIa. When measuring fu in a Monte-Carlo simulation, IIa would 
corresponds to the assumed underlying population, from which signals are randomly 
drawn in each test. The complementary detection probability rj = 1 — fu is 

/>oo 

ry(S*;II A )= / P df(6|H s H A ) d& . (12) 
Je, 

Note that this contains the usual definition of the power function t?(6*; A) as a special 
case, where IIa = A. The definition ( 12 1 can also be interpreted as the expected power 
over a population IIa- The Neyman-Pearson framework for hypothesis testing defines 
the most powerful test of size f a as a test that has the highest detection probability 
r] (i.e. smallest false dismissal fu) for a false-alarm probability of at most /a (e.g. see 
[23j[24]). In this framework one compares the respective detection probabilities rj of 
different detection statistics at a given false-alarm probability f^, which defines the 
receiver- operator characteristics (ROC), namely the function ^(/a; IIa). 

In general the relative performance of different detection statistics will depend on 
the assumed probability distribution IIa of signal parameters. One statistic can be 
more efficient in certain parts of parameter space and less efficient in others. If a test 
is most powerful over the whole parameter space (i.e. has the highest r)(&*; A) for all 
A S A), it a called a uniformly most powerful test. 



4. Frequentist maximum-likelihood approach: the .F-statistic 

When comparing two simple hypotheses, such as 7^n and Tls(A) for known A, the 
Neyman-Pearson Lemma states that the most powerful test (cf. section 3.3 1 is the 
likelihood-ratio C, defined as 



C(x; A) = exp 



pdf(x\H s (A)) 
pdf(x|H N ) 
ing {: 

(x\s) 



C(x; A) 

Assuming Gaussian noise and using (10 1, we explicitly obtain 

1 



(s\s) 



(13) 



(14) 



However, in the case of a composite hypothesis H.s with unknown amplitude 
parameters A, the orthodox frequentist framework does not generally provide a 
canonical detection statistic. Interestingly, one cannot even define a frequentist 
likelihood pdf(x|7is). A common but ad-hoc approach to dealing with composite 
hypotheses consists of using the maximum of the likelihood ratio C(x; A) over the 
parameter space A, i.e. define 



max£(i: ^4) 



(15) 



as a composite-hypothesis test: decide for 7ig if Cml(x) > and H.n otherwise. 
Using (14 1 and Q, the likelihood-ratio function can be written more explicitly as 



C(x\ A) = exp 



-A^M^A" 



where we defined 



Xf, = (x|/i M (A s )) , = (MAsOIMAg)) • 



(16) 
(17) 
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We see that C(x;A) is a Gaussian function in A^, so we can analytically maximize it 
to obtain 

£ M LW=e rW . with T{x)= l -x^M^ (18) 

with the inverse matrix ./VP" defined via A4^ a M. av — 5£. This defines the so-called 
^-statistic, which was first derived in this context in [5]. The statistic 2T{x) can 
be shown to be ^-distributed with four degrees of freedom, and a non-centrality 
parameter p 2 = (s|s), where p is called the (optimal) signal-to-noise ratio. The 
expectation value of 2T is ^[2^] = 4 + p 2 . 



5. Bayesian hypothesis testing 



The Bayesian hypothesis-testing framework follows uniquely from a straightforward 
application of the probability axioms (cf [25l [26j [lOj El])- For any question of 
interest one can (at least in principle) compute the probability of different hypotheses, 
optimally using the available information such as the observed data x(t), and all 
our prior information and assumptions, which we denote by T. Here we use the 
Bayesian approach to construct a classical detection statistic, in order to compare its 
performance to the frequentist ^"-statistic in the Neyman-Pearson framework. 

For any hypothesis Hi we can directly express the probability of Hi being true 
given the data x and our background assumptions /, namely 

This expression is known as Bayes' theorem, and it follows directly from the product 
rule of probabilities applied to pdi(Hi x\I). The term P(Hi\ I) is the prior probability 
for Hi- Contrary to the orthodox frequentist framework, the (marginal) likelihoocjf] 
pdf(x\Hi I) of observing data x given Hi is well-defined even for composite hypotheses. 
In order to compute pdf(x\Hs I), we simply use the product rule to write 

pdf(^ s /)= pd{(xlHsI) , (20) 

and invoking the normalization condition J pdf (A\xHs I) d A = 1, we obtain 

pdf(x\H s I)= f pdfiAxlHs^d^A 

7 (21) 
= / pdf(x\AHsI)pdf(A\H s I)d 4 A. 

J A 

If {Hi}™ =1 is a set of m mutually exclusive and exhaustive hypotheses, i.e. exactly one 
of them is true, then one obtains the normalization condition P{7~Li\ x I) — 1; 



which determines the denominator pdf(x|J) in (19 1. We do not need to make this 
assumption, however, as we can instead compute the relative probability of H§ with 
respect to H^, which is known as the (posterior) odds ratio Osn, namely 

P(H s \xI) pdi(x\H$I) P(H S \I) 
SNl l) ~ P(H N \xI) pdf(x\H N I) P(H N \I) ' 1 1 



X also known as the evidence 
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This expressions shows how the prior odds ratio P(7is\I) / 'P(Wn|-0 gets "updated" by 
the observation of x, namely by multiplication with the (marginal) likelihood-ratio 

pdi( x \n s i) f . 

BsN(x|/) = pdf(,|W N 7) ' (23) 

which is also known as the Bayes factor. Note that the prior odds ratio is a constant 
factor in Osn ; an d therefore plays no role in constructing a classical detection statistic 



&(x) (any monotonic function of &(x) has the same power). Using (211 and (131, we 
can write the Bayes factor (p3| explicitly as 



B SN (x\I)= [ C(x; A) pdf (A\H S I) d 4 A. (24) 

J A 



Note that while the jF-statistic ( 18 ) was obtained by maximizing the likelihood ratio 
£(x;A) over the "nuisance parameters" A^ , the Bayes factor B SN (x) consists of 
marginalizing C{x;A) with an amplitude prior pdi(A\Hs I). In order to uniquely 
specify the Bayes factor, we therefore need to determine the function pdf (A\HsI), 
which adequately describes our ignorance of the pulsar amplitude parameters A. 



5.1. Uniform priors in A^ -coordinates: rediscovering the T -statistic 



Considering the form ( 16 1 of £(x;A), a straightforward choice would be uniform priors 



in coordinates A 11 . We refer to this as the "canonical prior" II C , namely 

pdf(M'}l« s n o /) = { c ;^ l,< """ TO 

where ft,jf ax is the maximum amplitude we consider possible, and ho({A^ 1 }) is 
given by inversion of Eq. (j5j. The normalization constant C is determined by 
f A pdf(A\Hs I) d 4 A = 1. The actual choice of h™ ax is unimportant for the properties 
of Bs^(x) as a detection statistic, because for large /i™ ax 3> 1, the marginalization 
( pl| leads to a Gaussian integral, namely 

B SN (x\U c I) = C [ C(x;A)d 4 A 

Ja (26) 

= C{2^) 2 {c\ctM)- 1/2 e :F(x \ 
where det M. is the determinant of the matrix M.a V . We see that uniform amplitude 



priors in ^-coordinates lead us back to the jF-statistic (18). However, there is an 
additional antenna-pattern factor (det M)^ 1 ^ 2 , which depends on the sky-position n 
and the observation period. For a targeted search with a single known sky position, 
this is a constant factor which does not affect the power of the detection statistic, i.e. 
-Bsn(^|II c /) is equivalent to T{x). This weighting factor would play a role, however, 
when investigating searches over unknown sky position. A similar effect was first noted 
in the Bayesian analysis of burst detection statistics [15] . 



5.2. Physical priors in amplitude- space: introducing the B-statistic 

Despite our assumed "ignorance" about the amplitude parameters of the GW signal, 
we have made a number of model assumptions about the geometry of the emitting 
system (see section |2| . In a sense the physical model describing the emitting system 
singles out a preferred coordinate system in A, in which we should express our 
ignorance. We refer to the resulting prior as the "physical prior" n p h- 
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The initial phase <po is directly related to the rotation angle of the quadrupolar 
deformation with respect to the rotation axis f2 at some fixed reference time. The 
probability distribution for (f>o is therefore independent of {ho, cos t, ip}, and by 
rotational symmetry we can assign a uniform prior, i.e. 

pdf(0o|n ph /) = i 0o S [0,2tt). (27) 
Ztt 

In section [2] and Fig. [I] we have seen that cos i and ip determine the orientation of 
the rotation axis £1 with respect to the observer frame ^, — fi}. If we have no 
information about the orientation of the emitting system, then rotational symmetry 
dictates an isotropic probability distribution for f2. The surface element on the unit 
sphere of orientations can be expressed as d 2 S = \dcos l\ \dip\, and because cos i and 
i\) are independent degrees of freedom, their respective prior probabilities are 

pdf (cos i |n ph I) = - , cos l e [-1, 1] , (28) 



P df(v|n ph /) 



2 
2 



ip € [-tt/4, tt/4) 



(29) 



where we used the fact that gauge transformations {ip — > ip + 7r/2, 0o — > 0o + 
leave the observed signal Q unchanged, so i/j can always be brought into the range 
ip G [— 7r/4, 7r/4). Note that these priors are identical to those used previously in 
Bayesian parameter estimation O [TBI C7] and Monte-Carlo simulations for frequentist 
upper limits [71 ITS]. 

Contrary to the angle variables cos i,ip and (f>o, there is no unique natural choice 
of uninformed prior for the amplitude ho- One could derive a prior for ho from 
Eq. ([2]), if astrophysical priors for the deformation s, spin rate f2 and distance d are 
available. Other possibilities include a "maximum entropy" prior, or a Jeffreys prior. 
For simplicity, however, we simply chose a uniform prior, namely 

1 



pdf(^ |n ph 7) 



u max 

n 



Combining d27l 



- ( 30 1 we obtain an amplitude prior of the form 
pdf(/i ,cost» |n ph 7) = 2?r2 1 femax = C , h £ [0,h 



max] 
J 



(30) 



(31) 



which for simplicity of notation we refer to as the physical prior, while this qualifier 
can only be justified for the angle variables. Substituting the prior H p h in the Bayes 
factor (24), and assuming h™* ^> 1, we now obtain 

B sn (x\Il ph I) = C / dh a dcosi dip d(f> C(x; A) , (32) 
Jo J -l J-f Jo 

with the likelihood-ratio £(x;A) of Eq. (16). We use this Bayes factor as a new 



classical detection statistic B(x), which we refer to as the "S-statistic", namely 

B{x)=B su {x\U ph I). (33) 



5.3. Relation between amplitude priors H c and n p h 



In order to compare the physical amplitude prior n p h of Eq. (31 1 to the canonical 
prior II C of Eq. ( 25 1 , we use the coordinate transformation (I5J) relating A^ and A 1 . 



The Jacobian J of this transformation is found as 



J = 



det (m 



h 3 
4 V 



2 \3 
COS t 



(34) 
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Using the identity 

pdf {{A") | n c I) d 4 A = pdf ({A 1 } | n c I) d 4 A , (35) 

together with the relation d 4 A = J d 4 A between volume elements, we can translate Il c 
into physical coordinates {A 1 } = {/io, cos t, -0, 0o}, namely 

C 3 
pdf (/i ,cos l, ip, 0o | n c I) = — h\ (l — cos 2 t) , (36) 



which can be compared to the physical prior n p h in Eq. (31 1. We see that Il c agrees 
with Ilph in assigning uniform prior probabilities to </>o and ip, but the prior densities 
on cost and h are very different. 

The canonical prior information Il c , which is implicit in the ^-statistic (cf. 
section 5.1), is therefore found to be rather unphysical: a higher prior probability 



is assigned to stronger signals compared to weaker ones, and signals with near- 
linear polarization (cost ~ 0, corresponding to "edge-on" emitters) are given undue 
weight compared to signals with near-circular polarization (| cosi| ~ 1, corresponding 
to "face-on" emitters). This amounts to postulating a non- isotropic probability 
distribution for the orientation f2 of spinning GW sources, such that Cl favors 
orientations orthogonal to the line of sight h. 



6. Comparing detection efficiencies of T(x) and B(x) 



6.1. Estimating the ROC curves 



We use the classical Neyman-Pearson framework (cf. section 3.3 1 in order to compare 
the detection efficiency, or "power", of the JF-statistic (18 1 and the ^-statistic (32). 
The practical Monte-Carlo procedure for estimating the ROC curve rj(j for 
any detection statistic 6 is straightforward: first generate a large sample of Nmc 
random draws {6n} of the statistic & for the case of no signal, i.e. s = 0. From 



this distribution we can estimate the false-alarm probability (111 as a function of the 
threshold 6* 



namely 
/a(6*) 



>s» 



Nmc 



(37) 



where NgyQ^ is the number of Sn values found above the threshold 6*. Similarly, in 
the signal case 7is, we randomly draw signal parameters from the assumed population 
II& and generate corresponding random draws {©s} of the statistic 6. From this 
distribution we can estimate the detection probability (12 1, namely 



77(6*; n A ) 



(38) 



Inverting (37 1 to yield 6* (/a) and substituting this into (38 1, we obtain the ROC 
curve 77(/a;IIa). 



6.2. Parameters used in Monte-Carlo simulation 

In targeted searches the Doppler parameters A of the signal are known, and for 
simplicity of this example we fixed these parameters as: right ascension a = 2rad, 
declination 5 = — 0.5rad and a constant frequency without spindown. We assume the 
detector location to be LIGO Hanford, and an observation with GPS start time of 
to = 756950413 and duration of T = 25 hours. The resulting numerical components of 
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the antenna-pattern matrix M.^ v of Eq. (17 1 are found as 

M22 = Ma = §-B, and M 12 = M 3i = jf~ C, 
C = —0.0104, with all other components (approximately) zero 



with A = 0.154, B = 0.234 
These parameters 

are given for the sake of completeness, the qualitative conclusions do not depend 
on these choices. We used Nmc = 10 6 random draws for each distribution, and we 
estimate the errors on rj(fA\ n&) using a jackknife estimator (see [28]) with 100 subsets. 
The estimated la errors on the detection probability in the following ROC curves are 
always less than a(rj) < 0.004. 



6.3. Monte-Carlo results 

Because Hs is a composite hypothesis, the ROC curves depend on the choice of 
injected signal population Ha. In order to illustrate the dependency on the amplitude 
parameter space, we first consider two highly unphysical choices of signal populations, 
namely (i) 14a consisting of a single, linearly polarized signal with cos 1 = 0, iji = 
and (ii) a single, (nearly) circularly polarized signal with cos 6 = 0.99, ip = 0. In 
both cases we fixed the SNR of the signal to be p = 4. Note that the choice of </>o is 
irrelevant for both T and B. These two choices reflect universes in which all spinning 
GW sources happen to be (i) edge-on or (ii) face-on, without the observer having 
any knowledge about it. The results of these simulations are shown in Fig. [2] For 
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Figure 2. ROC curves r/(fj\;Ilj,J comparing ^-statistic JlSj, B-statistic (33} , 
and the perfect-match likelihood ratio £(_4 S )- The chosen signal populations 11a 
have fixed SNR of p = 4, and consist of (i) a single linearly polarized signal with 
cos l = 0, V = (left panel), and (ii) a (nearly) circularly polarized signal with 
cos i = 0.99, V = (right panel). 



comparison purposes we also plot the perfect-match likelihood ratio C(A S ), which is 
optimal by the Neyman-Pearson lemma for testing the simple hypothesis Hs(A s ) (cf. 
section 



3.3 1, but requires all signal parameters ^4 S to be exactly known beforehand. 



We see that the ^F-statistic is more powerful than the ^-statistic if the signal 
is linearly-polarized, while the ^-statistic dominates for (near-) circularly-polarized 
GWs. This is not surprising given that the implicit ^F-statistic prior II C is biased in 
favor of linear polarization (cf. section 5.3 1. The results in Fig. [2] show that neither 
T- nor 23-statistic is uniformly most powerful (cf. section |3.3[ ) over the amplitude 
parameter space A. Note that this does not imply that the JF-statistic is more 



Bayesian versus maximum-likelihood statistics 



11 



powerful if we know a given source to be (near-) linearly polarized. One would fold 
this knowledge into the prior II in the Bayes factor p4| , while there is no natural way 
in which this knowledge can be incorporated into the ^-statistic. The resulting Bayes 
factor would therefore be more powerful than the ^"-statistic. 

In the next step we look at more realistic situations in which the injected signals 
are drawn from a population IIa of randomly distributed cos i and ip, according to 
the physical prior n p h (cf. section 5.2 1, with (iii) a fixed SNR of p = 4 and (iv) a 
fixed amplitude of h = 10 \/S^. The results of these simulations are shown in Fig. [3] 
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Figure 3. ROC curves r)(fA\^K) comparing ^-statistic JTSJ, S-statistic | |33| l, 
and the perfect-match likelihood ratio C(A S )- The signal populations IIa consist 
of randomly distributed cost and ip, according to the physical prior ilp^ (cf. 
section |5.2| l, with (iii) a fixed SNR of p = 4 and (iv) a fixed amplitude of 
ho = 10V5^. 



We see that in these situations the B-statistic is consistently more powerful than the 
jF-statistic. This is not surprising, given that the amplitude prior II C that is implicit in 
the ^-statistic differs substantially from the injected "real-world" isotropic probability 
distribution IIa on the orientation of fl. The B-statistic prior n p h, on the other hand, 
is consistent with IIa by construction (cf. section 5.3 1. In fact, it can be argued [T^] 
that the Bayes factor with a signal prior that is consistent with the population injected 
in the Monte-Carlo simulation is by construction optimal in the sense of the highest 
expected detection probability at given false-alarm probability. One would therefore 
not expect any other detection statistic to outperform B(x) in the simulations shown 
in Fig. [3] 



7. Conclusions 



We have shown that the maximum-likelihood ^-statistic can be interpreted as a Bayes 
factor with a simple, but unphysical, amplitude prior (and an additional unphysical 
sky-position weighting) . Using a more physical prior based on an isotropic probability 
distribution for the unknown spin-axis orientation of emitting systems, we obtain a 
new detection statistic, referred to as B-statistic. Monte-Carlo simulations for signals 
with random (isotropic) spin-axis orientations show that the B-statistic more powerful 
(in terms of its expected detection probability) than the ^"-statistic. 

The ^-statistic is therefore not "optimal" in the classical sense. However, the 
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^"-statistic sensitivity appears to be quite comparable to the ^-statistic (see Fig. [3]) 
in the range of parameters considered, while being computationally more efficient 
(there are no integrations required), and fully characterized by a known simple 
distribution. Interpreting it as a Bayes factor clarifies its role as a statement about 
relative probabilities of hypotheses, and allows one to use the ^-statistic within 
a fully Bayesian framework. Such a choice would be based on the simplicity and 
computational efficiency of the ^-statistic, despite the fact that it is not "optimal" . 
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